Differential Gene Expression (DGE) refers to the process of identifying and quantifying changes in gene expression levels between different biological conditions, such as healthy vs diseased states, treated vs untreated samples, different developmental stages or distinct cell types. Understanding DGE is fundamental in fields like molecular biology, genetics, and biomedical research, as it provides insights into how genes contribute to various biological processes and phenotypes.
Tumor Necrosis Factor (TNF) is a pro-inflammatory cytokine, crucial for the regulation of the immune system’s healthy inflammatory reactions. However, when in overproduction or dysregulation, TNF can contribute to chronic inflammation implicated in various autoimmune and inflammatory diseases, hence the now widespread use of anti-TNF drugs to treat such conditions, also known as TNF inhibitors (Monaco et al. 2014). Despite the success of anti-TNF agents in treating chronic pathological inflammatory reactions, little is known about their impact on the affected tissues at the transcriptome level. As an endeavor to discover the effects of such drugs on the expression levels of different genes, Karagianni and colleagues applied four different anti-TNF drugs on an established mouse model of inflammatory polyarthritis and collected a large number of independent biological replicates from the synovial tissue of healthy, diseased and treated animals (Karagianni et al. 2019). The dataset that pertains to these experiments is used in this workflow, which is part of the study’s computational analyses for identifying and clustering differentially expressed genes. Expression levels were detected and quantified with DNA microarrays, where detected fluorescence indicates the expression of a specific gene against a reference sample.
2 Exploratory Data Analysis
Exploratory Data Analysis (EDA) is a critical process in data science and statistics that involves analyzing and summarizing the main characteristics of a dataset, often using visual methods. The goal of EDA is to uncover patterns, spot anomalies, test hypotheses, and check assumptions with the help of summary statistics and graphical representations (Morgenthaler 2009).
boxplot(data, horizontal = T, las =1, cex.axis =0.5)
Figure 1: Boxplot of data before quantile normalization
Code
# Show head of dataframekable(head(data)) |>kable_styling(bootstrap_options =c("striped")) |>scroll_box(width ="100%", height ="100%") |>kable_classic()
Exp1Wt_1
Exp1Wt_2
Exp1Wt_3
Exp2Wt_1
Exp2Wt_2
Exp2Wt_3
Exp2Wt_4
Exp3Wt_1
Exp3Wt_2
Exp3Wt_3
Exp1Tg_1
Exp1Tg_2
Exp1Tg_3
Exp2Tg_1
Exp2Tg_2
Exp2Tg_3
Exp2Tg_4
Exp3Tg_1
Exp3Tg_2
Exp3Tg_3
Exp4Tg_1
Exp4Tg_2
Exp4Tg_3
Exp1Rem3_1
Exp1Rem3_2
Exp1Rem3_3
Exp1Rem6_1
Exp1Rem6_2
Exp1Rem6_3
Exp2Rem6_1
Exp2Rem6_2
Exp2Rem6_3
Exp2Rem6_4
Exp3Rem6_1
Exp3Rem6_2
Exp3Rem6_3
Exp2Hum6_1
Exp2Hum6_2
Exp2Hum6_3
Exp2Hum6_4
Exp3Hum6_1
Exp3Hum6_2
Exp3Hum6_3
Exp4Hum6_1
Exp4Hum6_2
Exp4Hum6_3
Exp2Enb6_1
Exp2Enb6_2
Exp2Enb6_3
Exp2Enb6_4
Exp4Enb6_1
Exp4Enb6_2
Exp4Enb6_3
Exp4Enb6_4
Exp4Enb6_5
Exp4Enb6_6
Exp2Cim6_1
Exp2Cim6_2
Exp2Cim6_3
Exp2Cim6_4
Exp4Cim6_1
Exp4Cim6_2
Exp4Cim6_3
Exp4Cim6_4
Exp4Cim6_5
Exp4Cim6_6
A1bg
3.78911
3.90081
3.72526
4.23290
4.14451
4.17460
4.19904
4.37631
4.35622
4.35167
3.78454
4.07887
4.13286
4.05313
4.20157
4.18684
4.27081
4.50461
4.29096
4.36422
3.74793
3.74200
3.79036
3.59538
3.96438
3.79821
3.55839
3.72886
3.73221
4.00460
4.06560
4.10832
4.14768
4.45705
4.40942
4.29107
4.23619
4.11532
4.04981
4.05459
4.38727
4.36555
4.49019
3.55204
3.45422
3.75162
4.11680
4.11782
4.09583
4.05194
3.70226
3.75635
3.82224
3.39433
3.74648
3.48890
4.11928
4.14522
4.10418
4.10931
3.61480
3.37258
3.58425
3.88982
3.28356
4.04937
A1cf
4.21654
4.17014
4.05211
4.60213
4.62320
4.65975
4.57693
4.15842
4.21211
4.18009
3.98771
4.11709
4.10684
4.68407
4.63325
4.73126
4.67542
4.09641
3.97244
4.10268
3.30531
3.29409
3.06133
4.29460
4.22492
4.08466
4.08205
4.13011
4.09484
4.62164
4.66620
4.61722
4.55999
4.38385
4.06813
4.13370
4.53864
4.75150
4.64363
4.68842
4.03449
4.62612
4.24503
3.39752
3.35704
3.31550
4.61107
4.71722
4.65451
4.68732
3.34686
3.26266
3.04503
3.16430
3.13805
3.10227
4.65152
4.67835
4.67399
4.61035
2.90180
3.27259
3.02554
3.12053
3.08666
3.24268
A2ld1
5.32714
5.62827
5.66799
8.32861
8.27225
8.27447
8.42341
8.22626
8.39943
8.57076
5.23156
5.29504
5.28384
7.90894
7.85233
8.14223
7.79835
7.90913
8.08984
7.95827
6.55196
6.51255
6.43711
6.01249
5.46742
5.54192
5.56968
5.70214
5.46825
8.26541
8.30650
8.22540
8.09633
8.30025
7.86106
7.86827
8.19100
8.37507
8.20061
8.27231
8.39751
8.30546
8.49839
6.69627
6.52055
6.64086
8.00233
8.18752
8.25757
8.03628
6.43991
6.62674
6.80411
6.40195
6.77547
6.49056
8.09591
8.09534
8.08303
8.07122
6.51077
6.56018
6.40339
6.71310
6.35072
6.61334
A2m
5.28545
6.49319
5.19705
6.34696
5.68818
7.02162
6.30722
7.24319
6.80225
6.69186
4.71382
5.03411
5.41249
6.57654
6.10905
6.11674
5.67339
6.93047
6.70562
6.65830
5.29271
5.47738
6.79552
6.19873
5.40119
6.48264
5.78780
5.32115
5.44562
6.88086
5.34267
6.33201
6.39798
6.45636
6.43533
6.19721
6.10152
6.31918
5.95871
6.76260
6.76862
6.14994
5.71077
5.23411
5.38092
5.36692
5.98963
5.63186
5.84172
5.34126
5.51386
5.58334
5.51965
5.63152
5.64672
6.46401
5.42400
6.97807
6.39452
5.24073
4.64791
4.64806
6.10691
5.04173
4.90970
4.32605
A3galt2
4.59905
4.59698
4.52439
4.85959
5.07224
5.01061
5.07838
5.62244
5.79039
5.66514
4.73714
4.90568
5.18567
4.95670
4.65065
4.95141
4.96971
5.23615
5.38050
5.49452
4.72958
4.78397
4.66467
4.41278
4.82166
4.49756
4.54898
4.53977
4.85098
4.89872
4.96186
5.01156
5.02533
5.50526
5.51928
5.52995
5.12319
5.15686
5.02892
4.91145
5.58904
5.62234
5.64090
4.49644
4.98682
4.62059
4.94060
5.10253
4.92040
4.80770
4.63297
4.79730
4.73491
4.63326
4.78748
4.69899
4.96660
5.02830
4.96265
5.03972
4.70114
4.64965
4.50466
4.70696
4.55098
4.64297
A4galt
7.73303
7.65939
7.93984
8.55175
8.48415
8.47890
8.22892
8.72356
8.71276
8.58030
7.22117
7.00438
6.80762
8.65125
8.49504
8.61091
8.60534
8.90632
8.84868
8.88411
8.87908
8.68815
8.59815
7.98948
7.60469
7.47936
7.60934
7.44681
7.36165
8.70350
8.61929
8.80555
8.47461
8.76373
9.06743
9.36881
8.55062
8.73755
8.65360
8.69017
9.22244
9.12786
9.25526
8.76674
8.93546
8.86313
8.66064
8.60910
8.63829
8.86686
8.73614
8.96106
8.77789
8.77312
8.97563
8.70684
8.83717
8.88996
8.98343
8.75343
8.80162
8.76917
8.70345
8.80999
8.86269
8.62107
3 Check for missing values
If there are any missing values, we will need to decide how to handle them, probably by removing the respective genes.
Code
# Check genes_data for missing valuescolSums(is.na(data))
Data distribution refers to how data values are spread or dispersed across a range of possible values. Understanding the distribution of the data is a key part of Exploratory Data Analysis (EDA), as it gives insights into the central tendency, variability, and shape of the dataset. It can help in detecting patterns, outliers, and important characteristics of the data, and also informs which statistical techniques to use. Thus, a detailed assessment of data distribution is a foundational step in unraveling the complexities of any dataset during the EDA process.
Central Tendency
Mean (Average): The sum of all values divided by the number of values. It provides a sense of the center of the data.
Median: The middle value when data is ordered. It’s less affected by extreme values (outliers) and provides a robust measure of the center.
Mode: The most frequent value in the data. For categorical data, this is often the only measure of central tendency.
Spread (Dispersion)
Range: The difference between the maximum and minimum values.
Variance: The average of the squared differences from the mean. It quantifies how much the data points deviate from the mean.
Standard Deviation: The square root of the variance. It provides a more interpretable measure of spread in the same units as the data.
Interquartile Range (IQR): The range within the middle 50% of the data (between the first and third quartiles). It’s a robust measure of spread that isn’t affected by outliers.
Code
# Adjust the layout and marginspar(mfrow =c(8, 9), mar =c(1, 1, 2.5, 1))for (col incolnames(data)) {plot(density(data[[col]]), main = col,xlab = col, col ="#009AEF", lwd =2)}
Figure 2: Data distribution before quantile normalization
The primary goal of quantile normalization is to align the statistical distributions of all samples so that they are the same. This technique assumes that all samples should have the same overall distribution of values (e.g., gene expression levels), meaning that any differences between samples that aren’t due to biological variation are artifacts that need to be corrected.
Steps of Quantile Normalization
Sort the Data:
For each sample (e.g., each column in a matrix of gene expression data), the data is sorted / ranked in ascending order.
Average the Sorted Values:
Across all samples, for each rank, the values are averaged. This creates a “reference” distribution, where each rank corresponds to the average value of that rank across all samples.
Replace Original Values:
Each value in the original dataset is replaced by the corresponding value from the reference distribution.
Output the Normalized Data:
The result is a new dataset where each sample has the same distribution of values, which helps remove any technical variations between samples.
Code
# Convert dataframe to matrixdata =as.matrix(data)# Normalize datadata =normalize.quantiles(data , copy=TRUE)# Convert matrix to dataframedata =data.frame(data)# Add column names to dataframecolnames(data) = Sample# Add row names to dataframerownames(data) = Gene
Code
boxplot(data, horizontal = T, las =1, cex.axis =0.5)
Figure 3: Boxplot of data after quantile normalization
Code
# Adjust the layout and margins par(mfrow =c(8, 9), mar =c(1, 1, 2.5, 1))for (col incolnames(data)) {plot(density(data[[col]]), main = col,xlab = col, col ="#009AEF", lwd =2)}
Figure 4: Data distribution after quantile normalization
6 Dimension Reduction Algorithms
6.1 Uniform Manifold Approximation and Projection (UMAP)
UMAP (Uniform Manifold Approximation and Projection) is a dimensionality reduction algorithm that is widely used in data science and machine learning. It is designed to help visualize high-dimensional data in a lower-dimensional space, typically 2D or 3D, while preserving the global and local structure of the data as much as possible.
Key concepts of UMAP
Manifold Learning: UMAP is based on the idea that high-dimensional data often lies on a lower-dimensional manifold, a geometry concept. UMAP aims to learn this manifold and then project the data onto a lower-dimensional space that best represents the structure of the data.
Topology: UMAP uses topological techniques to understand the shape of the data. It constructs a graph that represents the relationships between data points based on their proximity. This graph is then used to project the data into a lower-dimensional space.
Optimization: The final embedding is obtained by optimizing a cost function that balances the preservation of local structures (how points are grouped) and global structures (the overall layout).
Applications
Visualization: UMAP is often used to visualize complex datasets, such as in biology (e.g., gene expression data), where it helps to identify patterns and clusters.
Preprocessing: It can be used as a preprocessing step before other machine learning algorithms, such as clustering or classification, to reduce dimensionality while maintaining the integrity of the data structure.
Data Exploration: UMAP is useful for exploring and understanding large datasets, making it easier to identify trends, outliers, and relationships within the data.
Comparison to other techniques
t-SNE: UMAP is similar to t-SNE, another popular dimensionality reduction technique. However, UMAP is generally faster and better at preserving both global and local structures, making it more effective for large datasets.
PCA (Principal Component Analysis): PCA is another common method for dimensionality reduction but focuses on preserving linear relationships, which may not capture the complexity of data as effectively as UMAP, especially in cases of non-linear structures.
Code
# Keep only WT and TG sampleswt_tg_df = data[, 1:23]# After dataframe transposition columns must represent geneswt_tg_df =t(wt_tg_df)
Code
# UMAP dimension reduction for wt and tg sampleswt_tg_df.umap <-umap(wt_tg_df, n_components=2, random_state=15)# Keep the numeric dimensionswt_tg_df.umap <- wt_tg_df.umap[["layout"]]# Create vector with groupsgroup =c(rep("A_Wt", 10), rep("B_Tg", 13))# Create final dataframe with dimensions and group for plottingwt_tg_df.umap <-cbind(wt_tg_df.umap, group)wt_tg_df.umap <-data.frame(wt_tg_df.umap)# Plot UMAP resultsggplotly(ggplot(wt_tg_df.umap, aes(x = V1, y = V2, color = group)) +geom_point() +labs(x ="UMAP1",y ="UMAP2",title ="UMAP plot",subtitle ="A UMAP Visualization of WT and TG samples") +theme(axis.text.x =element_blank(),axis.text.y =element_blank(),axis.ticks =element_blank() ))
Figure 5: UMAP plot
6.2 Principal Component Analysis (PCA)
Principal Component Analysis (PCA) is a statistical technique used to simplify complex datasets by reducing their dimensionality while preserving as much variance as possible. The main idea behind PCA is to transform the original data into a new coordinate system where the greatest variance in the data is captured in the first few dimensions (called principal components). Here’s how PCA works step-by-step:
Covariance Matrix Computation
Calculate the covariance matrix of the data. The covariance matrix is a square matrix that shows the covariance between different features. If there are \(n\) features, the covariance matrix will be of size \(n×n\). Covariance gives an idea of how much two features vary together.
Eigenvalue and Eigenvector Calculation
Compute the eigenvalues and eigenvectors of the covariance matrix.
Eigenvectors represent the directions of the new feature space, and eigenvalues tell us how much variance is along each of these directions.
The eigenvector with the highest eigenvalue is the first principal component, which captures the most variance in the data. The second principal component is orthogonal to the first and captures the second most variance, and so on.
Sort Eigenvectors
Sort the eigenvectors by their corresponding eigenvalues in descending order. This ordering helps in deciding which principal components to keep. The components corresponding to the largest eigenvalues capture the most significant structure in the data.
Dimensionality Reduction
Choose the top \(k\) eigenvectors (where \(k\) is the number of dimensions you want to reduce your data to) and project the original data onto these eigenvectors. This transformation gives you the data represented in terms of the principal components, reducing the dimensionality from \(n\) to \(k\).
The resulting \(k\)-dimensional space retains most of the variability in the original \(n\)-dimensional space.
Projection of data
Finally, the original data is transformed into the new feature space defined by the selected principal components. This projection is done by multiplying the original data matrix by the matrix of eigenvectors.
Applications
Data compression: Reducing the number of features while retaining most of the original information.
Noise reduction: By keeping the most significant components, noise in the data (which might be in less significant components) can be reduced.
Visualization: Reducing data to 2 or 3 dimensions makes it easier to visualize, especially in high-dimensional datasets.
Preprocessing: Used as a preprocessing step before applying machine learning algorithms to improve performance and reduce computational cost.
plot_grid(fviz_pca_ind(wt_tg_df.pca, repel =TRUE, # Avoid text overlappinghabillage = group,label ="none",axes =c(1, 2), # choose PCs to plotaddEllipses =TRUE,ellipse.level =0.95,title ="Biplot: PC1 vs PC2") +scale_color_manual(values =c('#33cc00','#009AEF95')) +scale_fill_manual(values =c('#33cc00','#009AEF95')),fviz_pca_ind(wt_tg_df.pca, repel =TRUE, # Avoid text overlappinghabillage = group,label ="none",axes =c(1, 3), # choose PCs to plotaddEllipses =TRUE,ellipse.level =0.95,title ="Biplot: PC1 vs PC3") +scale_color_manual(values =c('#33cc00','#009AEF95')) +scale_fill_manual(values =c('#33cc00','#009AEF95')),fviz_pca_ind(wt_tg_df.pca, repel =TRUE, # Avoid text overlappinghabillage = group,label ="none",axes =c(2, 3), # choose PCs to plotaddEllipses =TRUE,ellipse.level =0.95,title ="Biplot: PC2 vs PC3") +scale_color_manual(values =c('#33cc00','#009AEF95')) +scale_fill_manual(values =c('#33cc00','#009AEF95')),# Visualize eigenvalues/variancesfviz_screeplot(wt_tg_df.pca, addlabels =TRUE,title ="Principal Components Contribution",ylim =c(0, 65), barcolor ="#009AEF95", barfill ="#009AEF95"),# Contributions of features to PC1fviz_contrib(wt_tg_df.pca, choice ="var", axes =1, top =14, color ="#009AEF95", fill ="#009AEF95"),# Contributions of features to PC2fviz_contrib(wt_tg_df.pca, choice ="var", axes =2, top =14, color ="#009AEF95", fill ="#009AEF95"),labels =c("A", "B", "C", "D", "E", "F"))
Figure 6: Plot PCA results
Code
wt_tg_df.pca <-data.frame("PC1"= wt_tg_df.pca$x[,1], "PC2"= wt_tg_df.pca$x[,2], "group"= group)# Plot PCA results# insert dataframe [1] , variables [2]-[3] and color groyp [4]ggplotly(ggplot( wt_tg_df.pca, aes(x= PC1 , y= PC2 , color= group ))+# try "geom_point" or "geom_line"geom_point()+# try "ggtitle" or "ggname"ggtitle("Two First Components of PCA") +theme(axis.text.x =element_blank(),axis.text.y =element_blank(),axis.ticks =element_blank()))
Figure 7: Plot two first components of PCA
7 Statistical Analysis
7.1 Group treatments in dataframe
ANOVA, or Analysis of Variance, is a statistical method used to compare means among three or more groups to determine if there are any statistically significant differences between them. It extends the t-test, which is typically used to compare the means of two groups, to situations where more than two groups are involved. ANOVA allows for the simultaneous assessment of variations within and between groups, enabling the identification of genes with expression patterns that are significantly different across experimental conditions. It is particularly useful in experiments where multiple groups or conditions are being tested simultaneously.
Hypotheses in ANOVA
ANOVA tests the following hypotheses:
Null Hypothesis (H₀): All group means are equal. (No significant difference between groups)
Alternative Hypothesis (H₁): At least one group mean is different from the others. (There is a significant difference between groups)
Partitioning Variance
ANOVA works by partitioning the total variance observed in the data into two components:
Between-Group Variance: The variation due to the interaction between the different groups. This measures how much the group means differ from the overall mean.
Within-Group Variance (Error Variance): The variation within each group. This measures how much individual observations differ from their group mean.
Mathematically, the total variance is expressed as:
ANOVA calculates an F-ratio by comparing the between-group variance to the within-group variance:
\[ F = \frac{\text{Between-Group Variance}}{\text{Within-Group Variance}} \]
A large F-ratio indicates that the between-group variance is large relative to the within-group variance, suggesting that the group means are significantly different.
A small F-ratio suggests that any observed differences in means are more likely due to random chance.
ANOVA table
The results of ANOVA are typically summarized in an ANOVA table, which includes:
Sum of Squares (SS): Measures of variation for both between-group and within-group components.
Degrees of Freedom (df): The number of values that are free to vary for each component.
Mean Squares (MS): Calculated by dividing the sum of squares by the corresponding degrees of freedom.
F-Statistic (F): The ratio of mean squares between groups to mean squares within groups.
Assumptions of ANOVA
For the results of ANOVA to be valid, several assumptions need to be met:
Independence: The observations within each group and between groups should be independent of each other.
Normality: The data in each group should be approximately normally distributed.
Homogeneity of Variances: The variances among the groups should be approximately equal.
Code
# -------------------Apply ANOVA on all genes----------------------# Create Matrix by Excluding rownames and colnamesmatrixdata =as.matrix(data)# Create Groupsgroup =factor(c(rep("A_Wt", 10),rep("B_Tg", 13),rep("C_Proph_Ther_Rem", 3),rep("D_Ther_Rem", 10),rep("E_Ther_Hum", 10),rep("F_Ther_Enb", 10),rep("G_Ther_Cim", 10)))# Create empty dataframeanova_table =data.frame()# Recursive parse all genesfor( i in1:length( matrixdata[ , 1 ] ) ) {# Create dataframe for each gene df =data.frame("gene_expression"= matrixdata[ i , ], "group"= group)# Apply ANOVA for gene i gene_aov =aov( gene_expression ~ group , data = df)# Apply tukey's post-hoc test on ANOVA results tukey =TukeyHSD( gene_aov , conf.level =0.99)# Vector calling Tukey's values tukey_data =c(tukey$group["B_Tg-A_Wt", 1], tukey$group["B_Tg-A_Wt", 4], tukey$group["C_Proph_Ther_Rem-A_Wt", 1], tukey$group["C_Proph_Ther_Rem-A_Wt", 4], tukey$group["D_Ther_Rem-A_Wt", 1], tukey$group["D_Ther_Rem-A_Wt", 4], tukey$group["E_Ther_Hum-A_Wt", 1], tukey$group["E_Ther_Hum-A_Wt", 4], tukey$group["F_Ther_Enb-A_Wt", 1], tukey$group["F_Ther_Enb-A_Wt", 4], tukey$group["G_Ther_Cim-A_Wt", 1], tukey$group["G_Ther_Cim-A_Wt", 4], tukey$group["C_Proph_Ther_Rem-B_Tg", 1], tukey$group["C_Proph_Ther_Rem-B_Tg", 4], tukey$group["D_Ther_Rem-B_Tg", 1], tukey$group["D_Ther_Rem-B_Tg", 4], tukey$group["E_Ther_Hum-B_Tg", 1], tukey$group["E_Ther_Hum-B_Tg", 4], tukey$group["F_Ther_Enb-B_Tg", 1], tukey$group["F_Ther_Enb-B_Tg", 4], tukey$group["G_Ther_Cim-B_Tg", 1], tukey$group["G_Ther_Cim-B_Tg", 4])# Append Tukey's data to dataframe anova_table =rbind( anova_table , tukey_data)}colnames(anova_table) <-c("Wt_Tg_diff", "Wt_Tg_padj","Wt_Rem_P_diff", "Wt_Rem_P_padj", "Wt_Rem_diff", "Wt_Rem_padj", "Wt_Hum_diff", "Wt_Hum_padj", "Wt_Enb_diff", "Wt_Enb_padj", "Wt_Cim_diff", "Wt_Cim_padj","Tg_Rem_P_diff", "Tg_Rem_P_padj", "Tg_Rem_diff", "Tg_Rem_padj", "Tg_Hum_diff", "Tg_Hum_padj", "Tg_Enb_diff", "Tg_Enb_padj", "Tg_Cim_diff", "Tg_Cim_padj")# Add rownames with gene namesrownames(anova_table) = Gene
8 Volcano Plot
Code
# -----------------Volcano plot preparation---------------------# Create variables for upregulated/downregulated genes and genes with no observed change in expression levelsupWT =0downWT =0nochangeWT =0# Filter Differential Expressed GenesupWT =which(anova_table[ , 1 ] <-1.0& anova_table[ , 2 ] <0.05)downWT =which(anova_table[ , 1 ] >1.0& anova_table[ , 2 ] <0.05)nochangeWT =which(anova_table[ , 2 ] >0.05| (anova_table[ , 1 ] >-1.0& anova_table[ , 1 ] <1.0 ) )# Create vector to store states for each genestate <-vector(mode="character", length=length(anova_table[,1]))state[upWT] <-"up_WT"state[downWT] <-"down_WT"state[nochangeWT] <-"nochange_WT"# Identify names of genes differentially expressed between wt and tggenes_up_WT <-c(rownames(anova_table)[upWT])genes_down_WT <-c(rownames(anova_table)[downWT])# Union of DEGs between wt and tgdeg_wt_tg <-c(genes_up_WT, genes_down_WT)# Subset dataframe based on specific degsdeg_wt_tg_df <-subset(data , Gene %in% deg_wt_tg)## Dataframe for volcano plotvolcano_data <-data.frame("padj"= anova_table[,2], "DisWt"= anova_table[,1], state=state)
Code
ggplot(volcano_data , aes(x = DisWt , y =-log10(padj) , colour = state )) +geom_point() +labs(x ="mean(Difference)",y ="-log10(p-value)",title ="Volcano Plot",subtitle ="Differentially Expressed Genes (WT vs TG)") +# Insert line to show cutoffgeom_vline(xintercept =c( -1 , 1 ),linetype ="dashed",color ="black") +# insert line to show cutoffgeom_hline(yintercept =-log10(0.05),linetype ="dashed",color ="black")
Figure 8: Volcano plot of differentially expressed genes
8.1 Notes on Volcano Plot Interpretation
A volcano plot is a type of scatter plot commonly used in high-throughput data analysis, particularly in genomics, proteomics, and transcriptomics. It’s particularly useful for visualizing the results of differential expression analyses, where you’re comparing two conditions to identify genes, proteins, or other features that are significantly differentially expressed.
Structure
X-axis (mean(Difference)): This represents the magnitude of change between two conditions. A mean change of 0 indicates no change, positive values indicate upregulation (more abundant in the condition being tested), and negative values indicate downregulation (less abundant).
Y-axis (−log10(p-value)): This axis represents the statistical significance of the observed changes, often using a p-value from a statistical test. The higher the value on this axis, the more significant the change. Since it is a negative log scale, higher values represent smaller p-values.
Interpretation
Significance: Points located at the top of the plot have high statistical significance. The further a point is from the origin along the y-axis, the more significant the result.
Magnitude: Points farther to the right (positive mean change) indicate genes/proteins/features that are upregulated, while points farther to the left (negative mean change) indicate downregulation.
Significant Features: Typically, the most interesting features are those that are both statistically significant (high on the y-axis) and have a large magnitude of change (far left or far right on the x-axis). These appear as points in the upper left or upper right corners of the plot.
Thresholds: Horizontal and vertical lines are often added to represent thresholds for significance. Points outside these thresholds are often colored differently to highlight significant changes.
9 Uniform Manifold Approximation and Projection (UMAP) after identifying Differentially Expresssed Genes
Code
# Subset dataframe based on specific degsdeg_wt_tg_df =subset(data , Gene %in% deg_wt_tg)deg_wt_tg_df = deg_wt_tg_df[,1:23]# After dataframe transposition columns must represent genesdeg_wt_tg_df =t(deg_wt_tg_df)
Code
# UMAP dimension reduction for wt and tg samplesdeg_wt_tg_df.umap =umap(deg_wt_tg_df, n_components=2, random_state=15)# Keep the numeric dimensionsdeg_wt_tg_df.umap = deg_wt_tg_df.umap[["layout"]]# Create vector with groupsgroup =c(rep("A_Wt", 10), rep("B_Tg", 13) )# Create final dataframe with dimensions and group for plottingdeg_wt_tg_df.umap =cbind(deg_wt_tg_df.umap, group)deg_wt_tg_df.umap =data.frame(deg_wt_tg_df.umap)# Plot UMAP resultsggplotly(ggplot(deg_wt_tg_df.umap, aes(x = V1, y = V2, color = group)) +geom_point() +labs(x ="UMAP1", y ="UMAP2", title ="UMAP plot", subtitle ="A UMAP Visualization of WT and TG samples (DEGs subset)") +theme(axis.text.x =element_blank(),axis.text.y =element_blank(),axis.ticks =element_blank()) )
Figure 9: UMAP plot of differentially expressed genes
Code
# Group wt and tg as character and not factorgroup =c(rep("A_Wt", 10), rep("B_Tg", 13) )# Dimension reduction with PCA for wt and tg dataframedeg_wt_tg_df.pca =prcomp(deg_wt_tg_df , scale. =FALSE)deg_wt_tg_df.pca =data.frame("PC1"= deg_wt_tg_df.pca$x[,1] , "PC2"= deg_wt_tg_df.pca$x[,2] , "group"= group)# Plot PCA resultsggplotly(ggplot(deg_wt_tg_df.pca , aes(x=PC1,y=PC2,color=group))+geom_point()+labs(x ="PC1", y ="PC2", title ="PCA plot", subtitle ="A PCA Visualization of WT and TG samples (DEGs subset)") +theme(axis.text.x =element_blank(),axis.text.y =element_blank(),axis.ticks =element_blank()))
Figure 10: PCA plot of differentially expresssed genes
9.0.1 Identify differentially expressed genes between transgenic animals and at least one therapy
Code
# Volcano plot dataframe preparation for DEGs from TG vs therapiesupTHER =0downTHER =0nochangeTHER =0# Filter genes based on mean diff and p_value between TG and therapiesupTHER =which((anova_table[,13] <-1.0& anova_table[,14] <0.05) | (anova_table[,15] <-1.0& anova_table[,16] <0.05) | (anova_table[,17] <-1.0& anova_table[,18] <0.05) | (anova_table[,19] <-1.0& anova_table[,20] <0.05) | (anova_table[,21] <-1.0& anova_table[,22] <0.05) )downTHER =which((anova_table[,13] >1.0& anova_table[,14] <0.05) | (anova_table[,15] >1.0& anova_table[,16] <0.05) | (anova_table[,17] >1.0& anova_table[,18] <0.05) | (anova_table[,19] >1.0& anova_table[,20] <0.05) | (anova_table[,21] >1.0& anova_table[,22] <0.05) )nochangeTHER =which( ( (anova_table[,13] >-1.0& anova_table[,13] <1.0) | anova_table[,14] >0.05) | ( (anova_table[,15] >-1.0& anova_table[,15] <1.0) | anova_table[,16] >0.05) | ( (anova_table[,17] >-1.0& anova_table[,17] <1.0) | anova_table[,18] >0.05) | ( (anova_table[,19] >-1.0& anova_table[,19] <1.0) | anova_table[,20] >0.05) | ( (anova_table[,21] >-1.0& anova_table[,21] <1.0) | anova_table[,22] >0.05) )# Create vector to store states for each genestate =vector(mode ="character", length =length(anova_table[,1]))state[upTHER] ="up_THER"state[downTHER] ="down_THER"state[nochangeTHER] ="nochange_THER"# Identify names of genes differentially expressed between TG and therapiesgenes_up_THER =c(rownames(anova_table)[upTHER])genes_down_THER =c(rownames(anova_table)[downTHER])deg_tg_ther =c(genes_up_THER, genes_down_THER)# Combine DEGs from TG and THERDEGs =c(deg_tg_ther, deg_wt_tg)# Data frame with all DEGs for clusteringDEGsFrame = anova_table[rownames(anova_table) %in% DEGs, ]DEGsFrame =as.matrix(DEGsFrame)
# ---------------------Prepare data for heatmap-------------------------# Define custom function to perform hierarchical clustering with the Ward.D2 linkage methodhclustfunc =function(x)hclust(x, method ="ward.D2")# Define custom function to calculate pairwise Euclidean distances between data pointsdistfunc =function(x)dist(x, method ="euclidean")# Perform clustering on rows and columns cl.row =hclustfunc(distfunc(DEGsFrame[, c(1,3,5,7,9,11)]))# Extract cluster assignments of rowsgr.row =cutree(cl.row, k=6)# Apply a set of color palettecolors =brewer.pal(5, "Set3")heatmap <-heatmap.2( DEGsFrame[, c(1,3,5,7,9,11)],col =bluered(100), # blue-red color palettetracecol="black",density.info ="none",labCol =c("TG", "REM_P", "REM", "HUM", "ENB","CIM"),scale="none", labRow="", vline =0,mar=c(6,2),RowSideColors = colors[gr.row],hclustfun =function(x) hclust(x, method ='ward.D2'))
Figure 12: Heatmap of differentially expressed genes after hierarchical clustering
10.1 Notes on Heatmap Interpretation
Rows: Each row represents a gene and the genes are clustered hierarchically (dendrogram on the left). Genes with similar expression patterns are grouped together.
Columns: They represent experimental conditions/ sample types.
Color Gradient (Color Key):
Red represents upregulated genes (positive values).
Blue represents downregulated genes (negative values).
The scale ranges from -4 to +4, indicating the intensity of differential expression (e.g. mean change).
Clustering:
The dendrogram on the top shows how the conditions are related based on the gene expression patterns. Conditions that cluster closer together (e.g., REM and TG) are likely to have similar expression profiles.
The dendrogram on the left shows how genes are grouped based on similarity in expression across the conditions.
Key Differences of Hierarchical vs K-means Clustering
Aspect
Hierarchical Clustering
K-means Clustering
Approach
Builds hierarchy (Dendrogram)
Partitions into flat clusters
Number of Clusters
Determined post-analysis
Must be predefined
Cluster Shape
Can handle various shapes
Assumes spherical clusters
Distance Metric
Multiple metrics
(e.g. Euclidean, Manhattan)
Only Euclidean distance
Scalability
Not scalable for large datasets
Efficient and scalable
Visual Output
Dendrogram
Cluster assignment
Handling Outliers
Sensitive to outliers
Sensitive, but manageable
Performing k-means clustering followed by creating a heatmap with hierarchical clustering combines the strengths of both clustering methods to get a more insightful and comprehensive visualization of the data. This hybrid approach can be particularly useful when working with large and complex datasets, such as gene expression data. This combination leverages the strengths of both methods: k-means for global structure and hierarchical clustering for local, detailed exploration.
11 Functional Enrichment Analysis of Differentially Expressed Genes (DEGs)
Once DEGs are identified, they need to be annotated to determine their biological functions, cellular localization, molecular interactions, and involvement in various biological pathways. This is often achieved by comparing DEGs to databases of known gene annotations, such as Gene Ontology (GO) or Kyoto Encyclopedia of Genes and Genomes (KEGG).
Pathway analysis focuses on identifying interconnected networks of genes that collaborate to carry out specific biological functions or participate in common signaling pathways. This involves mapping DEGs onto existing biological pathways and identifying key regulatory nodes or hub genes within these pathways. Pathway analysis provides insights into the underlying molecular mechanisms driving the observed gene expression changes.
The code below performs hierarchical clustering analysis on a gene expression dataset and then further analyzes the clusters to identify enriched biological terms using the databases: Gene Ontology (GO), with three Sub-Ontologies (Biological Process (BP), Cellular Component (CC), Molecular Function (MF)) transcription factors (TF), and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases.
Code
# Get the cluster assignmentsmt <-as.hclust(heatmap$rowDendrogram)# Cut the tree into 8 clusterstgcluster <-cutree(mt, k =8)tgdegnames <-rownames(DEGsFrame)# Keep unique cluster numberscl <-as.numeric(names(table(tgcluster)))totalresults <-0totalcols <-0pcols<-c("firebrick4", "red", "dark orange", "gold","dark green", "dodgerblue", "blue", "magenta", "darkorchid4")# Iterating through clusters for functional enrichment # The gost() function queries the significant genes from each cluster against various biological databasesfor (i inc(6, 5, 4, 3, 2, 7, 1)) { gobp <-gost(query =as.character(tgdegnames[which(tgcluster == cl[i])]), organism ="mmusculus", significant = T, sources ="GO:BP")$result gomf <-gost(query =as.character(tgdegnames[which(tgcluster == cl[i])]), organism ="mmusculus", significant = T, sources ="GO:MF")$result gocc <-gost(query =as.character(tgdegnames[which(tgcluster == cl[i])]), organism ="mmusculus", significant = T, sources ="GO:CC")$result tf <-gost(query =as.character(tgdegnames[which(tgcluster == cl[i])]), organism ="mmusculus", significant = T, sources ="TF")$result kegg <-gost(query =as.character(tgdegnames[which(tgcluster == cl[i])]), organism ="mmusculus", significant = T, sources ="KEGG")$result# Combine results results <-rbind(kegg, tf, gobp, gomf, gocc)# Filter the results based on different sources tf <-grep("TF:", results$term_id) go <-grep("GO:", results$term_id) kegg<-grep("KEGG:", results$term_id)# Get enriched terms/pathways, their associated p-values, and other relevant information obtained from the enrichment analysis. kegg<- results[kegg, ] tf <- results[tf, ] go <- results[go, ]# Order the results based on p-values kegg<- kegg[order(kegg$p_value), ] go <- go[order(go$p_value), ] tf <- tf[order(tf$p_value), ]# Split the term_id and term_name ll <-strsplit(as.character(tf$term_name), ": ") ll <-sapply(ll, "[[", 2) ll <-strsplit(as.character(ll), ";") tf$term_name <-sapply(ll, "[[", 1)# Remove duplicatesif (length(tf$term_id) >0) { uniqtf <-unique(tf$term_name) tfout <-0for (ik in1:length(uniqtf)) { nn <-which(as.character(tf$term_name) ==as.character(uniqtf[ik])) tfn <- tf[nn, ] inn <-which(tfn$p_value ==min(tfn$p_value)) tfout <-rbind(tfout, head(tfn[inn, ], 1)) } tf <- tfout[2:length(tfout[, 1]), ] } results <-rbind(head(kegg, 10), head(go, 10), head(tf, 10)) totalresults <-rbind(totalresults, results) n <-length(results$term_id) totalcols <-c(totalcols, rep(pcols[i], n))}totalresults <- totalresults[2:length(totalresults[, 1]), ]totalcols <- totalcols[2:length(totalcols)]par(mar =c(5, 15, 1, 2))
Code
# Visualization of under-expressed clusters of DEGsbarplot(rev(-log10(totalresults$p_value[75:126])),xlab ="-log10(p-value)",ylab ="",cex.main =1.3,cex.lab =0.9,cex.axis =0.9,main ="Under-Expressed Clusters",col =rev(totalcols[75:126]),horiz = T,names =rev(totalresults$term_name[75:126]),las =1,cex.names =0.6)
Figure 13: Under-expressed clusters of differentially expressed genes after functional enrichment analysis
Karagianni, Niki, Ksanthi Kranidioti, Nikolaos Fikas, Maria Tsochatzidou, Panagiotis Chouvardas, Maria C. Denis, George Kollias, and Christoforos Nikolaou. 2019. “An Integrative Transcriptome Analysis Framework for Drug Efficacy and Similarity Reveals Drug-Specific Signatures of Anti-TNF Treatment in a Mouse Model of Inflammatory Polyarthritis.” Edited by Cinzia Cantacessi. PLOS Computational Biology 15 (5): e1006933. https://doi.org/10.1371/journal.pcbi.1006933.
Monaco, Claudia, Jagdeep Nanchahal, Peter Taylor, and Marc Feldmann. 2014. “Anti-TNF Therapy: Past, Present and Future.”International Immunology 27 (1): 55–62. https://doi.org/10.1093/intimm/dxu102.